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ABSTRACT 

The Phase Diverse Speckle (PDS) problem is formulated mathematically as Multi Frame Blind Deconvolution (MFBD) 
together with a set of Linear Equality Constraints (LECs) on the wavefront expansion parameters. This MFBD-LEC 
formulation is quite general and, in addition to PDS, it allows the same code to handle a variety of different data collection 
schemes specified as data, the LECs, rather than in the code. It also relieves us from having to derive new expressions 
for the gradient of the wavefront parameter vector for each type of data set. The idea is first presented with a simple 
formulation that accommodates Phase Diversity, Phase Diverse Speckle, and Shack-Hartmann wavefront sensing. Then 
various generalizations are discussed, that allows many other types of data sets to be handled. 

Background: Unless auxiliary information is used, the Blind Deconvolution problem for a single frame is not well 
posed because the object and PSF information in a data frame cannot be separated. There are different ways of bringing 
auxiliary information to bear on the problem. MFBD uses several frames which helps somewhat, because the solutions 
are constrained by a requirement that the object be the same, but is often not enough to get useful results without further 
constraints. One class of MFBD methods constrain the solutions by requiring that the PSFs correspond to wavefronts 
over a certain pupil geometry, expanded in a finite basis. This is an effective approach but there is still a problem of 
uniqueness in that different phases can give the same PSF. Phase Diversity and the more general PDS methods are special 
cases of this class of MFBD, where the observations are usually arranged so that in-focus data is collected together with 
intentionally defocused data, where information on the object is sacrificed for more information on the aberrations. The 
known differences and similarities between the phases are used to get better estimates. 

Keywords: Wavefront sensing, Deconvolution, Phase diversity, Inverse problems, Image restoration, Shack-Hartmann. 

1. INTRODUCTION 

We first present a technique for jointly estimating the common object and the aberrations in a series of images that differ 
only in the aberrations. With no extra information beyond the image formation model, including aberrations from the 
phase in the generalized pupil transmission function, it is a Maximum-Likelihood (ML) Multi-Frame Blind Deconvolution 
(MFBD) 1 method. Constraining the PSFs to be physical by requiring that they come from an underlying parameterization 
of the phase over the pupil is a powerful technique. Although such methods do work, 2 ' 3 methods using more information 
work better. 4 Data sets used in Phase-Diverse Speckle (PDS) interferometry 5 7 has such extra information in the form 
of two (or more) imaging channels with a known difference in phase (at least to type), with Phase Diversity (PD) 6 ' 8 ' 9 
as the special case of only one such pair (or set). These are all different data collection schemes lending themselves 
to similar joint inversion techniques. The purpose of the formulation presented here is to recognize the similarities and 
outlining a method for treating them all with a single algorithm. We discuss how the formulation accommodates also 
Shack-Hartmann (SH) wavefront sensing (WFS) as well as, with a couple of simple generalizations, several other types 
of data sets. 

We start with the ML metric for a MFBD data set and derive an algorithm for the joint estimation of object(s) and 
aberrations for a variety of imaging scenarios. In its simplest form, it is equivalent to MFBD, while the known relations 
between image channels for particular data collection schemes are expressed as Linear Equality Constraints (LECs). 

Throughout this paper a Gaussian noise model for the data is assumed because that permits simplifications in the ML 
solution methods. It is our experience with Gaussian noise PD, that it is a good model for low-contrast objects like the 
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solar photosphere 6 but also that artificial high-contrast objects (like a laser-pinhole source) can often be imaged in such a 
way that the noise is negligible. 10 

A code based on the formulation presented here has been successfully applied to PD, PDS, and MFBD data of various 
kind, both from the Swedish Vacuum Solar Telescope as well as several different simulated data sets. 



2. MULTI-FRAME BLIND DECONVOLUTION 
2.1. Forward Model and Error Metric 

We use an isoplanatic image formation model with Gaussian white noise, which means that we assume that the optical 
system can be characterized by a generalized pupil function, which can be written, for an image frame with number 
j e {1,...,/}, as 

P j =A j expU<p j } (1) 

where <f>j is the phase and Aj is a binary function that specifies the geometrical extent of the corresponding pupil. A data 
frame dj can then be expressed as the convolution of an object, /, and a point spread function (PSF), sj = | gf ~ l {Pj}\ 2 . In 
the Fourier domain we get 

Dj(u) = F(u) ■ Sj(u) + Nj(u), (2) 

where S j is the OTF, Nj is an additive noise term with Gaussian statistics and u is the 2-D spatial frequency coordinate. 
For brevity, we will drop this coordinate for the remainder of the paper. 

We parameterize the unknown phases by expanding them in a suitable basis, {if/ m }, allowing for a part of the phase, 6j, 
to be excepted from the expansion, 

M 

4>j = Oj + ^a jm iff m ; Vy. (3) 

m 

The Gaussian white noise assumption allows us to use the inverse Wiener filter estimate of the object, 
where 

j 

Q = 7ohj+J]\Sj\ 2 , (5) 

j 

to derive a metric in a form that does not explicitly involve the object 8 and that has been shown to correspond to a ML 
estimate of the phases. 9 With two regularization parameters, 1 1 this metric can be written as 



J IV-' n'C m 1 ./ 



\a jm \ 2 . (6) 



When minimizing L, the y D bj term in Q has the effect of establishing stability with respect to perturbations in the object. 
Its use in Eq. (4) suggests setting y bj to something proportional to the inverse of the signal to noise ratios of the image 
data frames. The other regularization parameter, y wf , stabilizes the wavefront estimates and can be set by examining the 
relation between L and the RMS of the wavefront. 12 For simplicity, this term is given under the assumption that the 
i// m are Karhunen-Loeve (KL) 13 functions, where X m is the expected variance of mode m. This should work reasonably 
well also for low-order Zernike polynomials. In the general case, the wavefront regularization term is a matrix operation 
involving the covariances. 11 

Note that, although presented in a PD setting, this metric has nothing to do with PD per se, it is just a MFBD metric. 
The PD part enters in the parameterization in earlier PD methods and, equivalently, in the LECs in this presentation. 
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2.2. Gradient and Hessian for Traditional Phase Diversity 

Efficient minimization of L requires the gradients and, for some methods, the Hessian with respect to the aberration 
parameters. We start by considering the traditional PD problem, for which the normal equations can be written as 

A PD -<to-b PD -0, (7) 

where the elements of a are the coefficients of an expansion of the single phase that is at the pupil of all diversity channels. 
The M elements of b PD can be expressed as the Euclidean inner products, {■ , •), of an expression for the gradient of L and 
the aberration basis functions taken over the definition area of the basis functions, 11 



C = 1-2% Im[p* %{ Pj KclTHF-Dj - \F\ 2 Sj}]}] , A + y wf ^. 



(8) 



An approximation of the M x M Hessian matrix A PD can be derived if Q is regarded as a fixed quantity. 11 The elements 
can then be written as 



Q 



where 



U u = 



A 

Q 



Di 



(9) 



(10) 



Other PD gradient and Hessian formulae are also available in the literature. 6 ' 9 We base our derivations on the ones 
given here, because they incorporate the regularization terms. 

2.3. Gradient and Hessian for MFBD 

The difference between MFBD and PD is that we in general do not know any relation between the wavefronts in the differ- 
ent channels. The PD gradient and Hessian expressions can then be simplified considerably by relaxing all dependencies 
between imaging channels, because we can then set all multiplications involving pupil quantities with different j to zero. 
This corresponds to J different pupils with independent phases. 

We also have to multiply with a factor V7 for each differentiating operation. This is because we are splitting each 
independent variable a m into J different aj m . The vector sum of J instances of identical dor ;m is V7 times larger than 
da m . Because this factor appears in the denominator of the derivatives, we have to correct the gradient and the Hessian by 
multiplying with and J, respectively. 

We lexicographically arrange all the or ;m in a single column vector, a, with N = JM elements, so that we can also 
refer to them as a n , with a single index n = (j - l)M + m, 



a = [an a l2 ... a JM ] = [ari . 
We can then write the normal equations of the MFBD problem as 

A MFBD . 6q , _ b MFBD ^ q 

where the inner product part of Eq. (8) simplifies to the individual terms of the summations, 



fc™ = (-2 V7lm[p* ZlpjKeirHFDj - \F\ 2 S ,}]}] , * m ) + y wf ^. 



(11) 



(12) 



(13) 



Due to the independence between channels, the N x N Hessian A MFBD is block-diagonal, 

[Ai 



.MFBD 



(14) 
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where each M x M matrix Xj has A^™ D as its element at (m', m). With some arithmetics, Eq. (9) can be simplified so 
that these elements can be written as 



4 mfbd _ dJ 



im[ YJ^} *{p, <r'{^v,,}} + p) *{ PJ rr'^v,,}}) + 



(15) 



where 



(16) 



3. LINEAR EQUALITY CONSTRAINTS 

3.1. Notation and Theory 

Now we add information about the data set in the form of LECs. The theory for solving optimization problems with LECs 
can be found in some text books on numerical methods, e.g. Ref. 14. The idea is that each constraint reduces the number 
of unknowns by one. We can transform a constrained optimization problem with N parameters and Nc constraints into an 
unconstrained problem with N - Nc unknowns. 

The constraints are given as a set of linear equations that have to be satisfied exactly, 

Cor-d = 0, (17) 

while minimizing L. With Nc linearly independent constraints, C is an Nc x N matrix, where Nc < N. This means 
Eq. (17) is under-determined and can therefore be solved exactly by several different a. All solutions to Eq. (17) can be 
written as 

a = a + Q 2 p, (18) 

where a is a particular solution to Eq. (17) and the N' = N - Nc column vectors of Q2 are an orthogonal basis of the null 
space of C. A particular solution can always be found by setting the N' last elements of a to zero and solving the reduced 
system 

C'a'-d = 0, (19) 

where C is the upper left Nc x Nc submatrix of C and a' is a vector containing the first Nc elements of a. Note, though, 
that all known differences can be incorporated into 6j, so we can set d = without loss of generality, and we always have 
the particular solution a = 0. We can find a basis for the null space by using an orthogonal decomposition of C, such as 
SVD or QR. With QR factorization, we write 

C T =Q R; Q = [Qi Q 2 ], (20) 

where Q2 is the rightmost N' columns of Q and T used as a superscript denotes matrix transpose. 

The constrained minimization problem in a can be transformed into an unconstrained minimization problem in a 
reduced set of variables, /?, an N' element vector of parameters, /3„>. The normal equations for the transformed problem 
are obtained by left-multiplying Eq. (12) with and substituting Q2#/? (Eq. (18), a = 0) for 6a, 

Q2 A MFBD Q2 • 8fi — Q2 • b MFBD — 0. (21) 

Once we have a solution for j8, we easily get the solution for a from Eq. (18). 

In order to efficiently express the LECs, we need to be able to distinguish between frames collected in different ways. 
We will therefore expand the j index to a set of two indices, k and f, so that j = k + (t - \)K e {1, . . . , TK). We will use 
k € {1, . . . , K) for simultaneous exposures in diversity channels and t € {1, . . . , 71 for discrete time or different realizations 
of atmospheric turbulence. With K — 1 we have MFBD, while T — 1 corresponds to PD. If K — T = 1, we have the BD 
problem. 
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Figure 1. Constraints and null space matrices for PD with only Figure 2. Constraints and null space matrices for PD with two 
known differences in phase. K = 4 and M = 12. White and unknown phase differences, e.g. the tilts. K = 4 and M = 12. 
black correspond to ± 1, respectively. See Eqs. (22) and (23). See Eq. (24). Compare Fig. 1. 



3.2. Phase Diversity with only Known Differences in Phase 



For a PD data set (T = 1) we know that the phases are equal, except for diversity and registration. The diversity is in 8k 
and does not enter in the expansion in basis functions. To reduce the KM wavefront parameters to M unknowns, we need 
a constraints matrix consisting of (K - \)M constraints, which can be written as 



-PD' 



aim - a km = 0; k > I, Vm, 



(22) 



where we momentarily disregard the registration of the imaging channels. The rows can be in any order but it seems 
natural to let m vary faster than k. Using Eq. (20), we get a null space matrix. These matrices can be written in block 
matrix form as 

Im ~Im flfl 



C PD ' = 



M 



I 



M 



M 



(23) 



where Im is the M x M identity matrix. C PD is a (K - l)M x KM-matiix and the corresponding null space matrix is 
KMxM, see example in Fig. 1. Note that this particular null space matrix (for every constraints matrix, there are infinitely 
many) is sparse, so matrix multiplications with Q2 can be fast, even for large problems. 

Q2 is easily interpreted by looking at Eq. (18). The K individual aj are identical copies of f)/ y[K. Identical aj are 
exactly what we expect, given the constraints. Comparing Eqs. (21) and (7), it is also easy to see that Ql ■ b MFBD = b PD 
(except for an inconsequential constant factor in the wavefront regularization term), so the same solution minimizes the 
problem in either formulation. For each specific K it is also easy to confirm that Q2A MFBD Q2 = A PD . So the entire 
original PD problem from Sect. 2.2 is retained. 

3.3. Phase Diversity with Partly Unknown Differences in Phase 

For modes with unknown differences between the channels, we simply do not add any constraint, thus allowing the 
algorithm to optimize the corresponding coefficients independently in the two channels based on the available data. This 
can be used e.g. for the focus mode, if the diversity it is not known well, or for other low order modes that may differ 
between diversity channels. 

The most common unknown differences, however, are for the tilt modes, corresponding to image registration. Unless 
we have information from other data sets, we generally don't know the tilt differences. So we want to exclude the tilt 
modes from Eq. (22). However, we may wish to prevent a common tilt to grow without bounds by requiring that they add 
up to zero, 

c pd . lljk akm = 0; me {tilt modes) 

I aim - atm — 0; k > 1, m £ {tilt modes). 



(24) 



This will still allow changes in the individual tilt coefficients. A sample ((K - 1)(M - 2) + 2) x KM constraints matrix 
and a corresponding null space matrix are shown in Fig. 2. The modes with unknown differences (tilts) are numbered as 
m e {1,2). The null space matrix generated by my QR factorization code is not very regular but it does work, however it 
is not as easily interpreted as the one in the previous subsection. Also, it is not as sparse. See discussion in Sect. 5. 
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Figure 3. Constraints and null space matrices for PDS with only Figure 4. Constraints and null space matrices for PDS with two 
known differences in phase. J = 4, K = 2 and M = 12. See unknown phase differences, e.g. the tilts. J = 4, K = 2 and 
Eq. (25), ignore special treatment of tilt modes. M = 12. See Eq. (25). 

The number of unknowns is the M - 2 non-tilt wavefront parameters for the common phase plus the registration of K 
channels, which amounts to 2(K-\) tilt parameters. This corresponds to the total number of MFBD wavefront parameters, 
KM, minus the number of constraints, (K - 1)(M - 2) + 2. 



3.4. Phase Diverse Speckle 

Again starting with the case of completely known phase differences, the number of wavefront parameters is KTM, but 
there are only T different atmospheric realizations. The real number of unknowns is therefore TM and we need (K- \)TM 
constraints. This is T times the number of constraints in the PD problem. For each t e {1, . . . , T] we simply add one set 
of constraints of the second type in Eq. (24), see sample matrices in Fig. 3. 

In practice, we must allow for unknown registration between PD channels as well as between consecutive PD sets. 
Normally we know that the difference in tilt between the PD channels, k, is constant over time, t. Finally, we require that 
the registration parameters sum up to zero independently in the two directions. 



-iPDS 



Ht,k a thn = 0; me {tilt modes) 

a t km ~ a t (k-\)m = a(t-i)km - a(t-i){k-\) m ; t>\, k > I, me {tilt modes) (25) 

octim - atkm = 0; W, k > 1, m $ {tilt modes). 

Note that this problem is just as constrained as the problem with completely known phase differences, thus the null space 
matrices have the same dimensions. See sample matrices in Fig. 4. The special tilt constraints result in a less orderly and 
sparse null space matrix, although not as bad as in the PD case. 

When working with seeing degraded data, it is likely that all PD sets are not of the same quality. One strength with 
PDS is that bad data can be helped by being processed jointly with the worse data. However, intuition suggests that there 
are cases where the bad data disturbs the inversion of the good data. It can then be useful to calculate individual metrics 
for each PD set, use Eq. (6) but skip the y W f term and limit the sum over j to one t. 



^Z[Z 



IA*I 2 - * ** 



(26) 



where the sum in Q, is limited correspondingly. This gives good diagnostics for determining which PD sets resulted in 
good wavefront estimates. 

3.5. Shack-Hartmann WFS 

A Shack-Hartmann (SH) wavefront sensor consists of a microlens array that form a number of sub-images, each sampling 
different parts of the pupil. SH wavefront sensing works with point sources as well as with extended objects. An example 
of the latter is that it has been used with solar telescopes. 15 ' 16 A fundamental limitation for such applications is that the 
number of microlenses is restricted by the requirement that each sub-image must resolve solar fine structure. 
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Usually, the mean gradient of the phase is estimated by measuring the relative or absolute image motion in the subim- 
ages. These local tilts are then combined into a phase over the whole pupil. This means that in conventional analysis 
of SH wavefront data, information about wavefront curvature over the microlenses, manifested as differential blurring in 
the subimages, is discarded. It is only by more detailed modeling that all relevant information on the wavefronts can be 
extracted, which in principle should be made down to the noise level of the data. 

With different Aj for different imaging channels, j, Eqs. (l)-(3) accommodate SH data. The phase over the full pupil 
is constrained to be the same for all channels, while the different Aj correspond to the different parts of the pupil that are 
sampled by the different microlenses. The constraints matrix is the same as for PD with J channels. 

Calibration of the relative positions of the SH sub-images is required for conventional processing. If such data are 
available, they can be entered as tilts in the different 6j. However, it should be possible to estimate these tilts from the data 
if 9j is set to zero and local tilts are included in the parameterization of the wavefront. On the other hand, no global tilts 
are necessary. 

It has been demonstrated in simulation with two microlenses that this image formation model allows higher modes 
to be estimated than with conventional local tilt approximations, because local aberrations within each subimage is taken 
into account. 17 SH data have also been considered for regular PD analysis. The microlens array acts as a single phase 
screen and the sub-images are regarded as a single phase diverse image, to be processed together with a conventional high- 
resolution image. 18 ' 19 Both approaches are sensitive to getting the geometry of the microlens array right (phase screen in 
the latter case, positions of the Aj in the MFBD-LEC case). However, in the MFBD-LEC approach, the registration of 
the subimages does not have to be known; it can be estimated. 

You can still use a conventional image as another image channel in the MFBD-LEC formulation. If local tilt calibra- 
tion data is available, this should be straight-forward. One would want to avoid estimating local microlens tilts for the 
high-resolution channel. This can easily be done by adding LECs that set them to zero. See also Sect. 4. 1 . 

4. EXTENDING THE PRESENTED FORMULATION 

The formulation presented so far was chosen so that it would be fairly straight-forward, while still accommodating PDS 
and SH WFS. In order to demonstrate the versatility of the LEC approach, we now discuss briefly a number of relaxations 
of some of the requirements, and the types of data sets that can then be treated. 

4.1. Different Phase Parameterization 

It is trivial to relax the requirements that M is the same for every data set. The single index of the a and b arrays can then 
be written as n = m + YJyl\ M f rather than n = JM. 

One application for this is the SH processing of Sect. 3.5. If the SH sub-images are to be processed together with a 
conventional high-resolution image, one may want to exclude the local sub-pupil tilts from the parameterization of the 
phase seen by the conventional camera. See also the following subsections. 

4.2. Different Objects 

We can also relax the requirement that all images have the same object. This requires involving another index, say 
s e {1, . . . , S } (for "scene" or "set"), and changing the formulae for the metric and its derivatives to allow for separate 
objects F s . The metric can then be written as L — L s , where L s is from Eq. (6) but summing over j only that involves 
object F s . 

One application for this is high-resolution solar magnetograms, see Ref. 20,21. Magnetograms can be made by 
calculating the difference between opposite polarization components of light in a Zeeman sensitive spectral line collected 
through a birefringent filter. Because the magnetic signal is essentially the difference between two very similar images, it is 
important to minimize artifacts from a number of error sources, among them registration, seeing, in particular differential 
seeing, and instrumental effects. 

The differential seeing problem can be handled by making each image the sum of many short exposures, where bad 
frames are discarded and the images can be well registered before adding. The short exposures can be made simulta- 
neously in the two polarization channels by using a polarizing beamsplitter. To reduce also the influence of different 
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aberrations in the paths after the beamsplitter, the polarization states can be rapidly switched between the detectors, so 
that each polarization state is recorded with a aberrations from both paths. 

In order to reduce the effects of seeing degradation of the image resolution, the two sequences of short exposures could 
also be MFBD restored for seeing effects. Better would be to do this using constraints that come from the fact that the 
aberration differences between the two channels is constant, with proper respect paid to the switching. Since registration 
parameters are estimated along with the phase curvature terms, this could solve the seeing problem while at the same 
time registering the image channels to subpixel accuracy. This should be better than registration with cross-correlation 
techniques, since non-symmetric instantaneous PSFs would tend to corrupt such alignment. 

The fact that the differences between the two polarization signals are so subtle might inspire another approach. If the 
object is considered to be the same for wavefront sensing purposes, then we can treat the data as a PDS data set with zero 
diversity terms 0j. The object is the same and so are the phases, except for the registration terms. 

4.3. Different wavelengths 

If PDS data of the same object is collected in several different wavelengths, and the cameras are synchronized, it has been 
shown that wavefronts estimated from one PDS set can be used for restoration of the object seen in the data collected 
in the other wavelengths. 22 However, this requires information about the wavefront differences as seen by the different 
cameras. The differences that come from the different wavelengths can easily be calculated and other differences can be 
estimated from a selected subset of a larger data set and then applied to all the data under the assumption that they do not 
change. 

This assumption could also be used to advantage in a joint treatment of the simultaneous PDS (SPDS) sets. The 
different-object generalization from Sect. 4.2 is used together with PDS LECs and LECs that express this assumption, 



where A s is the wavelength (in arbitrary units) used for set s. For the assumption of wavelength independent aberrations 
to hold, the optical system should be approximately achromatic. In practice, this means the differences between the 
wavelengths should not be too large. 

Again, this can also be relevant for SH, if a high-resolution image is processed together with the SH data - and the 
wavelength is not the same in the SH sensor as in the imaging camera. 

4.4. Different Numbers of Diversity Channels 

Finally, we can also treat data sets with different numbers of phase diversity channels in different wavelengths. This 
requires a K that varies with s, J — Y, s K s . This is perhaps most useful because it permits setting K s = 1 for some s. This 
corresponds to PDS data set in one wavelength and simultaneous MFBD data in another wavelength as in Refs. 23, 24. 

We should be able to run several such sets jointly and constrain the aberration differences between different s to be the 
same for all f, without requiring the object to be the same. 

Again, the advantage of the joint treatment is of course that wavefront differences between different imaging cameras 
do not have to be calibrated but are estimated together with the aberrations that are in common to the different cameras. 



We have presented a formulation of the MFBD problem that accomodates the varying data collection schemes involved in 
PD and PDS as well as SH WFS and a number of combinations of data in different wavelengths or polarizations and with 
or without diversity in the phase or pupil geometry. In doing this, we have not exhausted the types of data sets that can be 
treated with the method presented in this paper. As long as there are known relations between the imaging channels, that 
can be formulated in terms of constraints on the wavefront coefficients, the data could be handled and the constraints can 
be used to benefit wavefront sensing and image restoration. 



Tjs,t,k a stkm = 0; 

a stkm ~ a st(k-l)m = a ' s(t-\)km ~ a s(t-\)(k-\)m\ 
Ostlm ~ Wstkm = 0; 



V*, m e {tilt modes) 
V*, t > 1, k > 1, me {tilt modes) 
V*, Vf, k > 1, mi {tilt modes) 
s > 1 , t > 1, mi {tilt modes) 



(27) 



5. DISCUSSION 
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We have written a code that implements the method presented in this paper. It has been tested for PD, PDS and MFBD, 
but the other suggested processing strategies have not been tried yet. Trying some of the suggested strategies will be the 
subject of future papers. 

The code is new and as we apply it to different problems we anticipate significant improvements in the implementation 
of the method as well as in the method itself. One particular area of improvement is the calculation of the null space matrix. 
It consists of basis vectors for the null space and for any constraints matrix, there are infinitely many null space matrices, 
corresponding to rotations of the coordinate system in the null space. A method for making a null space matrix that is as 
sparse and regular as possible would be very useful. Matrix multiplications involving sparse matrices can be performed 
much faster than full matrices when the size of the problem is large. Also, a null space matrix such as the one in Fig. 1 is 
so much more instructive and easy to interpret than the one in Fig. 2. 
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